clear all
% This file simulates changes in the policies to incentivice the purchas
% eof EV.
%%%%%%%%%%%%%%%%%%%%%%%%%%% This version considers constant total sales
%%%%%%%%%%%%%%%%%%%%%%%%%%% across simulations to address only substitution
%%%%%%%%%%%%%%%%%%%%%%%%%%% between inside options
%first version 07/15/2021
%last version 10/04/2021
load('results_final')


theta1=theta1_updated;
theta2=theta2_updated;
ob=length(p);
xi=xi_updated;
omega=omega_updated;
mc=mc_update;

% Calculate expected purged cost
mc_0=mc_update./(1+tariffs);
delta_rec=x_d*theta1_updated(1:k_D);

num_sim=100;

% Create supply matrix characteristics switching production location


policies={'none','trade tariff','tax exemption','both'}; % policy can be 'tax exemption' 'trade tariff' 'none'
results={};

%%
for jj=1:length(policies)
    
    policy=policies{jj};
    new_sales_tax=ones(size(p))*0.19;
    new_tariffs=ones(size(p))*0.35;

    if strcmp(policy,'none')
        cf_salestax=sales_tax;
        cf_tariffs=tariffs;
    elseif strcmp(policy,'tax exemption')
        cf_salestax=new_sales_tax;
        cf_tariffs=tariffs;
    elseif strcmp(policy,'trade tariff')
        cf_salestax=sales_tax;
        cf_tariffs=new_tariffs;
    elseif strcmp(policy,'both')
        cf_salestax=new_sales_tax;
        cf_tariffs=new_tariffs;
    else
        print('enter valid policy')
    end

% set new price vector
    p_cf=(p_nt.*(1+cf_salestax))./1000;
    mc_bl2=mc_0.*(1+cf_tariffs);
    
    tic
    [sj_out, p_out] = new_eq(p_cf,mc_bl2,theta2,ns,ind_market,x2,cf_salestax, delta, A1, QI,year);
    toc
    results{jj}=[sj_out, p_out];
end

%% Calcualte the expected baseline variable profit
%save('output_sim_4_24.mat','results','-v6') % Last time saved 04/24
%load('sim_res.mat')
%load 'output_sim.mat'
%
none_res=results{1};
tariff_res=results{2};
tax_res=results{3};
both_res=results{4};

sim_res_ms=[none_res(:,1),tariff_res(:,1),tax_res(:,1),both_res(:,1)];
sim_res_p=[none_res(:,2),tariff_res(:,2),tax_res(:,2),both_res(:,2)];

v_c=getranddraw(length(theta2)-1, ns);
expmu1=expmu(theta2,v_c,x2,(p_nt.*(1+new_sales_tax))./1000);
sj_ob=mktshares(delta,expmu1,QI,ind_market);

sim_res_ms(year<2017,:)=repmat(sj_ob(year<2017),1,4);

%% Matching fiscal expenditure to no tariff excemption policy

tax_ava=find((sales_tax==0.05));% find vehicles with 5% sale tax
s_tar=0.35;
p_cf=(p_nt.*(1+sales_tax))./1000; % Baseline prices
mc_cf=mc_0; % Base line marginal cost
%s_mkt=sim_res_ms(:,1);

sum_ms=accumarray(ind_market,sj); % Total market share inside option
total_units=accumarray(ind_market,q); % Total units sold in each month
total_market=total_units./sum_ms;
sim_units=round(sim_res_ms.*total_market(ind_market),0);
[ww, zz] = ndgrid(ind_market,1:size(sim_units,2));

month_uni=accumarray([ww(:) zz(:)],sim_units(:));
month_ms=accumarray([ww(:) zz(:)],sim_res_ms(:));

total_market2=round(month_uni(:,1)./month_ms,0);

sim_units2=round(sim_res_ms.*total_market2(ind_market,:),0);

ST=sim_units2.*(sim_res_p-(sim_res_p./[1+sales_tax,1+sales_tax,1+new_sales_tax,1+new_sales_tax]));
TT=sim_units2.*(mc_0.*[tariffs,new_tariffs,tariffs,new_tariffs]);
FC=ST+TT;
%p_ini=(sim_res_p(:,2)./(1+sales_tax));
p_ini=sim_res_p(:,1)./(1+sales_tax);
ss=2;

[xx, yy] = ndgrid(year-2014,1:size(FC,2));

year_ST=accumarray([xx(:) yy(:)],ST(:));
year_TT=accumarray([xx(:) yy(:)],TT(:));
year_FC=accumarray([xx(:) yy(:)],FC(:));

[ww, zz] = ndgrid(ind_market,1:size(FC,2));

month_ST=accumarray([ww(:) zz(:)],ST(:));
month_TT=accumarray([ww(:) zz(:)],TT(:));
month_FC=accumarray([ww(:) zz(:)],FC(:));


save('CL_4_24.mat','p_ini','mc_cf','theta2','ns','ind_market','x2','tax_ava','s_tar', 'delta', 'A1', 'QI','year','FC','total_market2','ss','sj_ob')%Last time saved on april 2024
%%
[sj_out2, p_out2, tar_s] = new_eq_sj(p_ini,mc_update,theta2,ns,ind_market,x2,tax_ava,s_tar, delta, A1, QI,year,FC,total_market,ss);



expmu2 = expmu(theta2,v_c,x2,p_out2);
sj_out2=mktshares(delta,expmu2,QI,ind_market);

sj_out2(year<2019,:)=sj_ob(year<2019);
%save('output.mat')
%% Carbon Tax
tier=zeros(size(p_nt));
UVT=9.37;

tier(p_nt<=1363*UVT,:)=0.015;
tier((1363*UVT<p_nt)&(p_nt<=3066*UVT),:)=0.025;
tier(p_nt>3066*UVT,:)=0.035;

tier(id_fuelcat=='Electric',:)=0.01;
tier(id_fuelcat=='Hybrid',:)=0.01;

qd=[0.02:0.02:1];

quan=quantile(emissions,qd);

[temp,cat] = histc(emissions,quan);
cat(cat==50,:)=49;
cat=(cat+1).*1.5* UVT;

tic
[sj_out3, p_out3] = new_eq_CT(p_nt,mc_update,theta2,ns,ind_market,x2,sales_tax, delta, A1, QI,tier,cat,year);
toc

sj_out3(year<2019,:)=sj_ob(year<2019);
save results_sim2_04_24
%% Carbon Tax Second Scenario
p_ini2=sim_res_p(:,1)./(1+sales_tax);

tier2=zeros(size(p_nt));
UVT=9.37;
ctt=1.5; % Multiplication parameter for emissions categories.

sc=2; % Scenario for comparation

qd=[0.02:0.02:1];

quan=quantile(emissions,qd);

[temp,cat2] = histc(emissions,quan);
cat2(cat2==50,:)=49;
cat2=(cat2+1);
%save('CLCT.mat','sim_res_p','sales_tax','p_nt','emissions','mc_update','theta2','ns','ind_market','x2', 'delta', 'A1', 'QI','year','FC','total_market2','tariffs')

tic
[sj_out4, p_out4,ctt_c] = new_eq_CT_2(p_ini2,mc_update,theta2,ns,ind_market,x2,sales_tax, delta, A1, QI,year,FC,total_market2,tier2,cat2,UVT,sc,tariffs);
toc


%% plots variation in MS

sim_res_ms=[sim_res_ms,sj_out2,sj_out3,sj_out4];
sim_res_p=[sim_res_p, p_out2, p_out3/1000,p_out4];

sim_res_ms(year<2017,:)=repmat(sj_ob(year<2017),1,7);

tot_ms_sim=[];
fuel_cat=[];
s0_2=[];
mrkt_2=[];
mean_pr=[];

for qq=1:length(policies)+3
    temp=[];
    temp2=[];
    if qq==1
        for jj=1:max(ind_market)
            ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),sim_res_ms(ind_market==jj,qq));
            tempp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),sim_res_p(ind_market==jj,qq),[],@mean);
            tot_ms_sim=[tot_ms_sim;ttpp]; 
            fuel_cat=[fuel_cat;unique(id_fuelcat(ind_market==jj))];
            s0_2=[s0_2;repmat(1-unique(s0(ind_market==jj)),size(ttpp))];
            mrkt_2=[mrkt_2;repmat(jj,size(ttpp))];
            mean_pr=[mean_pr;tempp];
        end
    else
        for jj=1:max(ind_market)
            ttpp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),sim_res_ms(ind_market==jj,qq));
            tempp=accumarray(grp2idx(id_fuelcat(ind_market==jj)),sim_res_p(ind_market==jj,qq),[],@mean);
            temp=[temp;ttpp];
            temp2=[temp2;tempp];
        end
        tot_ms_sim(:,qq)=temp;
        mean_pr(:,qq)=temp2;
    end
end


day=ones(size(year));
date_mrk=unique(datetime(year,month,day));

plot(date_mrk,tot_ms_sim([1:4:288],1)./s0_2([1:4:288]))

hold on
plot(date_mrk,tot_ms_sim([1:4:288],2)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([1:4:288],3)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([1:4:288],4)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([1:4:288],5)./s0_2([1:4:288]))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','fix fiscal cost','Location', 'southoutside')
title('Market Share Variation in Electric Cars for Different Policies')

plot(date_mrk,tot_ms_sim([4:4:288],1)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([4:4:288],2)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([4:4:288],3)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([4:4:288],4)./s0_2([1:4:288]))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Market Share Variation in Internal Combustion Cars for Different Policies')

plot(date_mrk,tot_ms_sim([3:4:288],1)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([3:4:288],2)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([3:4:288],3)./s0_2([1:4:288]))
hold on
plot(date_mrk,tot_ms_sim([3:4:288],4)./s0_2([1:4:288]))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Market Share Variation in Hybrid Cars for Different Policies')

 %% Plot variation price
plot(date_mrk,mean_pr([1:4:288],1))
hold on
plot(date_mrk,mean_pr([1:4:288],2))
hold on
plot(date_mrk,mean_pr([1:4:288],3))
hold on
plot(date_mrk,mean_pr([1:4:288],4))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Price Variation in Electric Cars for Different Policies (1000 USD)')

plot(date_mrk,mean_pr([4:4:288],1))
hold on
plot(date_mrk,mean_pr([4:4:288],2))
hold on
plot(date_mrk,mean_pr([4:4:288],3))
hold on
plot(date_mrk,mean_pr([4:4:288],4))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Price Variation in Internal Combustion Cars for Different Policies (1000 USD)')

plot(date_mrk,mean_pr([3:4:288],1))
hold on
plot(date_mrk,mean_pr([3:4:288],2))
hold on
plot(date_mrk,mean_pr([3:4:288],3))
hold on
plot(date_mrk,mean_pr([3:4:288],4))
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','Location', 'southoutside')
title('Price Variation in Hybrid Cars for Different Policies (1000 USD)')

%% Consumer surplus

CS_none = CS(theta2,v,x2,sim_res_p(:,1),delta,ind_market);
CS_tariff = CS(theta2,v,x2,sim_res_p(:,2),delta,ind_market);
CS_tax = CS(theta2,v,x2,sim_res_p(:,3),delta,ind_market);
CS_both = CS(theta2,v,x2,sim_res_p(:,4),delta,ind_market);
CS_fcf = CS(theta2,v,x2,sim_res_p(:,5),delta,ind_market);
CS_ct = CS(theta2,v,x2,sim_res_p(:,6),delta,ind_market);
CS_ct2 = CS(theta2,v,x2,sim_res_p(:,7),delta,ind_market);

CS_table=[CS_none, CS_tariff, CS_tax,CS_both,CS_fcf,CS_ct,CS_ct2]*1000;

plot(CS_none)
hold on
plot(CS_tariff)
hold on
plot(CS_tax)
hold on
plot(CS_both)
hold on
plot(CS_fcf)
hold off;
legend('status quo','no trade tariff','no tax excemption','no policy','constant FC','Location', 'southoutside')
title("Consumers' Surplus Variation")

save results_sim5_24_04
%% 

plot(date_mrk,tar_s)
hold on
plot(xlim,[nanmean(tar_s(date_mrk>='01-Jan-2019')) nanmean(tar_s(date_mrk>='01-Jan-2019'))], 'r')
legend('sale tax','mean sale tax','Location', 'southoutside')
title("Counterfactual sale tax to match observed market share")
hyp_tar=tariffs;
hyp_tar(hyp_tar<0.35)=nanmean(tar_s(date_mrk>='01-Jan-2019'));

% p_cf2=(p_out2.*(1+new_sales_tax));
% mc_bl3=mc_0.*(1+hyp_tar);
%     
% tic
% [sj_out22, p_out22] = new_eq(p_cf2,mc_bl3,theta2,ns,ind_market,x2,new_sales_tax, delta, A1, QI,year);
% toc
% 
% sim_res_ms(:,5)=sj_out22;
% sim_res_p(:,5)=p_out22;

%% Tables

sim_units=round(sim_res_ms.*total_market(ind_market),0);
[ww, zz] = ndgrid(ind_market,1:size(sim_units,2));
[xx, yy] = ndgrid(year-2014,1:size(sim_units,2));

year_ms=accumarray([xx(:) yy(:)],sim_res_ms(:));
year_uni=accumarray([xx(:) yy(:)],sim_units(:));
month_uni=accumarray([ww(:) zz(:)],sim_units(:));
month_ms=accumarray([ww(:) zz(:)],sim_res_ms(:));

total_market2=round(month_uni(:,1)./month_ms,0);

sim_units2=round(sim_res_ms.*total_market2(ind_market,:),0);



%tot_mrk_tab=kron(total_market,ones(4,7)); % Market size table (4 scenarios, 4 fueltypes)
% tot_num_sim=round(tot_ms_sim.*tot_market2,0); % Table number of vehicles by scenario and by fueltype
% 
year_table=kron(unique(year),ones(48,1));
fuel_type_tab=kron(ones(72,1),[1:4]');

num_year_fc=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Num_unit_year=accumarray([year-2014 grp2idx(id_fuelcat)],sim_units2(:,qq));
    num_year_fc(:,qq,1)=Num_unit_year(end-1,:)';
    num_year_fc(:,qq,2)=Num_unit_year(end,:)';
end


Emi_unit_year=accumarray([year-2014 grp2idx(id_fuelcat)],emissions,[],@mean); % Average emissions by year and fuel cat


price_year_fc=zeros(4,4,2); %table with num of vehicles per year
for qq=1:6
    P_unit_year=accumarray([year_table fuel_type_tab],mean_pr(:,qq),[],@mean);
    price_year_fc(:,qq,1)=P_unit_year(end-1,:)';
    price_year_fc(:,qq,2)=P_unit_year(end,:)';
end

% Yearly market share sales

scenarios_ms=num_year_fc./sum(num_year_fc);% ms per scenario
 
sales=q.*p;

yearly_q=accumarray([year-2014,grp2idx(id_fuelcat)],sales);
yearly_tq=repmat(accumarray(year-2014,sales),1,4);

yearly_ms=yearly_q./yearly_tq;

%subplot(2,1)
bar(unique(year),yearly_ms(:,[1,3]))
legend('Electric','Hybrid', 'Location', 'southoutside')
xlabel('Year')
ylabel('Market sahre')

% Yearly q sold
yearly_q=accumarray([year-2014,grp2idx(id_fuelcat)],q);
yearly_tq=repmat(accumarray(year-2014,q),1,4);

yearly_q=yearly_q./yearly_tq;
clrs = [0 0 0; 0.9 0.9 0.9];

subplot(2,1,1)
char1=bar(unique(year),yearly_ms(:,[1,3]))
xlabel('Year')
ylabel('Market share')
set(char1,{'FaceColor'},{clrs(1,:),clrs(2,:)}.')

subplot(2,1,2); 
char2=bar(unique(year),yearly_q(:,[1,3]))
xlabel('Year')
ylabel('Share of units sold')
set(char2,{'FaceColor'},{clrs(1,:),clrs(2,:)}.')

hL = legend([char1],{'Electric','Hybrid'});
% Programatically move the Legend
newPosition = [0.9 0.49 0.1 0.05];
newUnits = 'normalized';
set(hL,'Position', newPosition,'Units', newUnits);


% Mean price
yearly_p=accumarray([year-2014,grp2idx(id_fuelcat)],p,[],@mean);
yearly_p(yearly_p==0)=mean(yearly_p([1,3],4));

plot(unique(year),yearly_p(:,[1,3,4]))
xticks(unique(year))
m = {'+','o','*','x','s'};
legend('Electric','Hybrid','Gasoline', 'Location', 'southoutside')
xlabel('Year')
ylabel('MSRP + taxes (1000 USD)')

hyp_tar=new_tariffs; % hypothetic tariff for fixed fiscal cost counterfactual
sel_mrkt=ind_market(year>=2017);
for ii=min(sel_mrkt):max(sel_mrkt) %these loop allocates tariffs to vehicles types
    temp=find(ind_market==ii);
    tar_1=hyp_tar(temp);
    [~,tax_ava_i]=ismember(tax_ava,temp);
    tax_ava_i(tax_ava_i==0)=[];
    tar_1(tax_ava_i)=tar_s(ii);
    hyp_tar(temp)=tar_1;
end

%hyp_tar(hyp_tar<0.35)=nanmean(tar_s(date_mrk>='01-Jan-2019'));

cat_t2=cat2.*ctt_c(ind_market)* UVT;


sim_units=round(sim_res_ms.*total_market2(ind_market,:),0);
sim_units(:,5)=sim_res_ms(:,5).*total_market2(ind_market,2);
sim_units(:,7)=sim_res_ms(:,7).*total_market2(ind_market,2);

sale_tax_mat=(sim_res_p-(sim_res_p./[1+sales_tax,1+sales_tax,1+new_sales_tax,1+new_sales_tax,1+new_sales_tax,1+sales_tax,1+sales_tax]));
sale_tax_mat(:,7)=(((sim_res_p(:,7)-(cat_t2./1000))./(1+sales_tax)).*(sales_tax));
sale_tax_mat(:,6)=(((sim_res_p(:,6)-(cat./1000))./(1+sales_tax+tier)).*(sales_tax+tier));
variation_saletax=sale_tax_mat-sale_tax_mat(:,1);

[xx, yy] = ndgrid(year-2014,1:size(variation_saletax,2));
income_saletax=accumarray([xx(:) yy(:)],variation_saletax(:));
tariff_mat=(mc_0.*[tariffs,new_tariffs,tariffs,new_tariffs,hyp_tar,tariffs,tariffs]);

other_tax=zeros(size(tariff_mat));
other_tax(:,6)=((((p_out3-(cat./1000))./(1+sales_tax+tier)).*tier.*year>=2019)+(cat./1000));% Tax as stated in Colombian legislation proposal


other_tax(:,7)=cat_t2./1000;% Counterfactual carbon tax coefficient (assuming no addicional tax on prices)


variation_tariffs=tariff_mat-tariff_mat(:,1);
income_tariffs=accumarray([xx(:) yy(:)],variation_tariffs(:));

fiscal_cost=sim_units.*(tariff_mat+sale_tax_mat+other_tax)*1000;
fiscal_cost_tab=accumarray([xx(:) yy(:)],fiscal_cost(:));
fiscal_cost_dif=(fiscal_cost_tab-fiscal_cost_tab(:,4));


% Change in emissions
ave_km=20000; % assuming a total driving of 20k km per year

%sim_units=round(sim_res_ms.*total_market2(ind_market,:),0);
tot_emi=(emissions.*sim_units2.*ave_km)/(10^6);
life_emi=tot_emi*10; % 10 years Lifetime emissions in metric tons;
variation_emi=life_emi-life_emi(:,4);
yearly_emi_var=accumarray([xx(:) yy(:)],variation_emi(:)); % Emissions in metric tons

yearly_emi=accumarray([xx(:) yy(:)],life_emi(:)); % Emissions in metric tons
yearly_emi_dif=yearly_emi-yearly_emi(:,4);



cost_emi=fiscal_cost_dif./yearly_emi_dif; % cost in USD por ton of CO2
cost_emi(isnan(cost_emi))=0;
var_cost_emi=cost_emi-cost_emi(:,1);

% Change in energy consumption
%life_elec=tot_elec*10;
%yearly_elec=accumarray([xx(:) yy(:)],life_elec(:)); % Emissions in metric tons
%yearly_elec_dif=yearly_elec-yearly_elec(:,4);

% Additional emissions from electricity consumption
%yearly_add_em=yearly_elec/(10^6)*211.09;
%yearly_add_mei_dif=yearly_add_em-yearly_add_em(:,4);
%yearly_emi_add=yearly_emi+yearly_add_mei_dif;
%yearly_emi_dif_add=yearly_emi_add-yearly_emi_add(:,4);

% Profit per vehicle
p_notax=(sim_res_p./[1+sales_tax,1+sales_tax,1+new_sales_tax,1+new_sales_tax,1+new_sales_tax,1+sales_tax,1+sales_tax]);
p_notax(:,7)=((sim_res_p(:,7)-(cat_t2./1000))./(1+sales_tax));
p_notax(:,6)=((sim_res_p(:,6)-(cat./1000))./(1+sales_tax+tier));

profit_v=(p_notax-(mc_0.*[1+tariffs,1+new_tariffs,1+tariffs,1+new_tariffs,1+hyp_tar,1+tariffs,1+tariffs]))*1000;

sim_units2=round(sim_res_ms.*total_market(ind_market,:),0);
profit_tot=sim_units.*profit_v;
profit_tot2=sim_units2.*profit_v;
profit_avg=accumarray([xx(:) yy(:)],profit_v(:),[],@mean); % Profit per vehicle

year_uni2=accumarray([xx(:) yy(:)],sim_units(:));
var_profit_avg=profit_avg-profit_avg(:,4);
PR_yearly=accumarray([xx(:) yy(:)],profit_tot(:)); % Total Profit per vehicle


PR_yearly2=accumarray([xx(:) yy(:)],profit_tot2(:)); % Total Profit per vehicle
var_PR_yearly=PR_yearly-PR_yearly(:,4);
var_PR_yearly2=PR_yearly2-PR_yearly2(:,4);
tot_units=accumarray([xx(:) yy(:)],sim_units2(:));

%% Welfare
year_month=accumarray(ind_market,year,[],@mean);
[zz, vv] = ndgrid(year_month-2014,1:size(CS_table,2));

CS_yearly=accumarray([zz(:) vv(:)],CS_table(:),[],@mean); % yearly consumer surplus
var_CS_yearly=CS_yearly-CS_yearly(:,4);
Total_CS=(CS_yearly.*tot_units);
vsr_TS=Total_CS-Total_CS(:,4);
%CS_per_con=CS_yearly./tot_units;

total_wel=PR_yearly+(CS_yearly.*tot_units); % yearly welfare
total_wel_2=total_wel+fiscal_cost_dif;
total_wel_var=total_wel-total_wel(:,4);
total_wel_var_2=total_wel_2-total_wel_2(:,4); % Used in tables
%var_total_wel=total_wel-total_wel(:,1); % Variation in total welfare



wel_emi=total_wel_var_2./yearly_emi_dif;

% weigthed average emissions per scenario


temp_emi=sim_units.*emissions;


fuel_type_tab=kron(ones(72,1),[1:4]');

emi_year_fc=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Emi_unit_year=accumarray([year-2014 grp2idx(id_fuelcat)],temp_emi(:,qq));
    emi_year_fc(:,qq,1)=Emi_unit_year(end-1,:)';
    emi_year_fc(:,qq,2)=Emi_unit_year(end,:)';
end

num_year_fc2=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Num_unit_year2=accumarray([year-2014 grp2idx(id_fuelcat)],sim_units(:,qq));
    num_year_fc2(:,qq,1)=Num_unit_year2(end-1,:)';
    num_year_fc2(:,qq,2)=Num_unit_year2(end,:)';
end

emi_year_type=emi_year_fc./num_year_fc2; % Weigthed average of emissions by vehicle type

tot_marlk_e=repmat(sum(num_year_fc2),4,1,1);

ms_year_emissions=num_year_fc2./tot_marlk_e;



% Market share and emissions by price
qd2=[0.2:0.2:1];

quan2=quantile(p,qd2);

[temp,p_cat] = histc(p,quan2);

p_cat=p_cat+1;

pr_year_fc=zeros(6,7,2); 
for qq=1:7
    pr_unit_year=accumarray([year-2014 p_cat],temp_emi(:,qq));
    pr_year_fc(:,qq,1)=pr_unit_year(end-1,:)';
    pr_year_fc(:,qq,2)=pr_unit_year(end,:)';
end

num_year_fc3=zeros(6,7,2); 
for qq=1:7
    Num_unit_year3=accumarray([year-2014 p_cat],sim_units(:,qq));
    num_year_fc3(:,qq,1)=Num_unit_year3(end-1,:)';
    num_year_fc3(:,qq,2)=Num_unit_year3(end,:)';
end

emi_year_p=pr_year_fc(1:end-1,:,:)./num_year_fc3(1:end-1,:,:); % Weigthed average emissions by price quintile


% Total emissions per year dist by quintiles

em_year_fc=zeros(6,7,2); 
for qq=1:7
    pr_unit_year=accumarray([year-2014 p_cat],tot_emi(:,qq));
    em_year_fc(:,qq,1)=pr_unit_year(end-1,:)';
    em_year_fc(:,qq,2)=pr_unit_year(end,:)';
end

tot_marlk_p=repmat(sum(num_year_fc3),6,1,1);

ms_year_price=num_year_fc3./tot_marlk_p;

num_year_sg=zeros(8,7,2); 
for qq=1:7
    Num_unit_year_sg=accumarray([year-2014 grp2idx(id_staseg)],sim_units(:,qq));
    num_year_sg(:,qq,1)=Num_unit_year_sg(end-1,:)';
    num_year_sg(:,qq,2)=Num_unit_year_sg(end,:)';
end

emi_year_sg=zeros(8,7,2); 
for qq=1:7
    Num_emi_year_sg=accumarray([year-2014 grp2idx(id_staseg)],temp_emi(:,qq));
    emi_year_sg(:,qq,1)=Num_emi_year_sg(end-1,:)';
    emi_year_sg(:,qq,2)=Num_emi_year_sg(end,:)';
end
emi_year_sg=emi_year_sg./num_year_sg;

tot_marlk_sg=repmat(sum(num_year_sg),8,1,1);

ms_year_sg=num_year_sg./tot_marlk_sg;

% Mean provit

temp_pi=sim_units.*profit_v;

num_year_fc4=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Num_unit_year4=accumarray([year-2014 grp2idx(id_fuelcat)],temp_pi(:,qq));
    num_year_fc4(:,qq,1)=Num_unit_year4(end-1,:)';
    num_year_fc4(:,qq,2)=Num_unit_year4(end,:)';
end

profit_year_type=num_year_fc4./num_year_fc2; % Weigthed average of emissions by vehicle type

%% Individual consumer surplus
wel_per_vehi=profit_tot+CS_table+fiscal_cost;

CS_none_ind = CS_ind(theta2,v,x2,sim_res_p(:,1),delta,ind_market);
CS_tariff_ind = CS_ind(theta2,v,x2,sim_res_p(:,2),delta,ind_market);
CS_tax_ind = CS_ind(theta2,v,x2,sim_res_p(:,3),delta,ind_market);
CS_both_ind = CS_ind(theta2,v,x2,sim_res_p(:,4),delta,ind_market);
CS_fcf_ind = CS_ind(theta2,v,x2,sim_res_p(:,5),delta,ind_market);
CS_ct_ind = CS_ind(theta2,v,x2,sim_res_p(:,6),delta,ind_market);
CS_ct2_ind = CS_ind(theta2,v,x2,sim_res_p(:,7),delta,ind_market);

CS_table_ind=[CS_none_ind, CS_tariff_ind, CS_tax_ind,CS_both_ind,CS_fcf_ind,CS_ct_ind,CS_ct2_ind]*1000;


temp_cs=sim_units.*CS_table_ind;

cs_year=zeros(4,7,2); %table with num of vehicles per year
for qq=1:7
    Num_cs_year=accumarray([year-2014 grp2idx(id_fuelcat)],temp_cs(:,qq));
    cs_year(:,qq,1)=Num_cs_year(end-1,:)';
    cs_year(:,qq,2)=Num_cs_year(end,:)';
end

cs_year_type=cs_year./num_year_fc2; % Weigthed average of emissions by vehicle type

num_year_fc3=zeros(6,7,2); %table with num of vehicles per year
for qq=1:7
    Num_unit_year3=accumarray([year-2014 p_cat],sim_units(:,qq));
    num_year_fc3(:,qq,1)=Num_unit_year3(end-1,:)';
    num_year_fc3(:,qq,2)=Num_unit_year3(end,:)';
end

cs_year_fc=zeros(6,7,2); 
for qq=1:7
    cs_unit_year=accumarray([year-2014 p_cat],temp_cs(:,qq));
    cs_year_fc(:,qq,1)=cs_unit_year(end-1,:)';
    cs_year_fc(:,qq,2)=cs_unit_year(end,:)';
end

cs_year_type_q=cs_year_fc./num_year_fc3; % Weigthed average of emissions by vehicle type

ele_year_fc=zeros(6,7,2); % Number of electric vehicles per Q
for qq=1:7
    ele_unit_year=accumarray([year-2014 p_cat],(sim_units(:,qq).*(id_fuelcat=='Electric')));
    ele_year_fc(:,qq,1)=ele_unit_year(end-1,:)';
    ele_year_fc(:,qq,2)=ele_unit_year(end,:)';
end

tot_year_fc=zeros(6,7,2); 
for qq=1:7
    tot_unit_year=accumarray([year-2014 p_cat],sim_units(:,qq));
    tot_year_fc(:,qq,1)=tot_unit_year(end-1,:)';
    tot_year_fc(:,qq,2)=tot_unit_year(end,:)';
end
ms_ele_q=ele_year_fc./tot_year_fc;

hyb_year_fc=zeros(6,7,2); % Number of hybrid vehicles per Q
for qq=1:7
    hyb_unit_year=accumarray([year-2014 p_cat],(sim_units(:,qq).*(id_fuelcat=='Hybrid')));
    hyb_year_fc(:,qq,1)=hyb_unit_year(end-1,:)';
    hyb_year_fc(:,qq,2)=hyb_unit_year(end,:)';
end

ms_hyb_q=hyb_year_fc./tot_year_fc;

gas_year_fc=zeros(6,7,2); % Number of gasoline vehicles per Q
for qq=1:7
    gas_unit_year=accumarray([year-2014 p_cat],(sim_units(:,qq).*(id_fuelcat=='NR')));
    gas_year_fc(:,qq,1)=gas_unit_year(end-1,:)';
    gas_year_fc(:,qq,2)=gas_unit_year(end,:)';
end

ms_gas_q=gas_year_fc./tot_year_fc;

temp_p=sim_res_p.*sim_units;

gas_p_fc=zeros(6,7,2); % Prices of gasoline vehicles per Q
for qq=1:7
    gas_unit_year2=accumarray([year-2014 p_cat],(temp_p(:,qq).*(id_fuelcat=='NR')));
    gas_p_fc(:,qq,1)=gas_unit_year2(end-1,:)';
    gas_p_fc(:,qq,2)=gas_unit_year2(end,:)';
end

p_gas_q=gas_p_fc./tot_year_fc;

% Welfare by vehicle

wel_per_vehi=profit_tot+temp_cs+fiscal_cost;

wel_emi_vehi=wel_per_vehi;%./life_emi;
%wel_emi_vehi(isinf(wel_emi_vehi))=0;
%p_cat2=p_cat;
%p_cat2(isinf(wel_emi_vehi))=[];

wel_year=zeros(6,7,2); %table with cost_emi per quintile
for qq=1:7
    Num_unit_year3=accumarray([year-2014 p_cat],wel_emi_vehi(:,qq),[],@mean);
    wel_year(:,qq,1)=Num_unit_year3(end-1,:)';
    wel_year(:,qq,2)=Num_unit_year3(end,:)';
end

emi_year=zeros(6,7,2); %table with cost_emi per quintile
for qq=1:7
    Num_unit_year3=accumarray([year-2014 p_cat],life_emi(:,qq),[],@mean);
    emi_year(:,qq,1)=Num_unit_year3(end-1,:)';
    emi_year(:,qq,2)=Num_unit_year3(end,:)';
end

wel_emi_q=wel_year./emi_year;